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ABSTRACT 

Aims. The hlack hole candidates exhibit fast variability of their X-ray emission on a wide range of timescales. If the variability is 
recurrent, but not strictly periodic, we can find excess of power about certain frequency, often with Lorentzian shape. These peaks are 
referred to as quasi-periodic oscillations and they range between mHz up to kHz. Their origin is still under debate and lots of different 
scenarios has been proposed, including resonant effects in the black hole accretion flow. The purely stochastic variability that occurs 
due to turbulent conditions in the plasma, is quantified by the power density spectra and appears practically in all types of sources and 
their spectral states. The specific kind of quasi-periodic flares is expected, when the global structure of the accretion flow, governed 
by the nonlinear hydrodynamics, induces fluctuations around a fixed point solution. These conditions, which occur at high accretion 
rates, lead to the variability in the sense of deterministic chaos. One of plausible mechanisms proposed to explain such variability is 
the intrinsic thermal-viscous instability in the accretion disk. 

Methods. We study the nonlinear behaviour of X-ray sources using the recurrence analysis method. We estimate quantitatively the 
indications for deterministic chaos, such as the Renyi’s entropy, for the observed time series, and we compare them with the surrogate 
data. This powerful method, widely known in other fields of physics, is used for the first time in the astrophysical context. 

Results. Using the data collected by RXTE satellite, we reveal the oscillations pattern and the observable properties of six black hole 
systems. For five of them, we confirm the signatures of deterministic chaos being the driver of their observed variability. 
Conclusions. Both the well known microquasar GRS 1915-1-105, as well as its recently discovered analogue, IGR 117091-3624, 
exhibit variability characteristic to deterministic chaotic system. Therefore the underlying nature of the process must be intrinsically 
connected in these sources with the accretion flow instability, that leads to the limit cycle oscillations around a fixed point. Further¬ 
more, we found significant traces of non-linear dynamics also in three other sources (GX 339-4, XTE J1550-564 and GRO11655-40), 
particularly in the disk dominated soft state, as well as in the intermediate states at the rising and declining phase of the outburst. 
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1. Introduction 

The high energy radiation emitted by black hole X-ray binaries 
originates in an accretion disk, where the viscous dissipation of 
the gravitational potential energy is responsible for both heat¬ 
ing of the plasma and transporting of the angular momentum. 
Hardly any of the observed black hole systems exhibits the radi¬ 
ation flux being constant with time. Instead, most of the sources 
undergo fast and complicated variability patterns on different 
timescales. The variations that are purely stochastic in their 
nature, are expected since the viscosity of the accretion disk is 
connected with its turbulent behaviour induced by magnetic in¬ 
stabilities. Furthermore, when the variations occur close to some 
frequency, the excess in the power spectrum known as the quasi- 
periodic oscillations (QPOs) appears. The mechanism of the 
origin of QPOs is generally not known and different scenarios 
of their origin has been proposed, including oscillations of ac¬ 
cretion disc/torus ( |Titarchuk & Osherovich|2000||Chakrabarti &| 
|Manickam 2000|l or of shock f ronts inside the accreting material 
( |Das112003 [Sukova & Janiuk1|2015| l and non-linear resonances 
between the radial and orbital epicyclic motions ( |Abramowicz 
|et al.|[2006| ), possibly also relating the high frequency and low 
frequency QPOs (see the review by [Belloni & Stella| ( |2014) l). 
Among the non-resonant class of models, the Lense-Thirring 
precession is often invoked and aimed to explain the twin kHz 


QPOs as well as the low frequency QPOs, and was originally 
proposed for the neutron star X-ray binaries in their Z and Atoll 
classes ( [Stella & Vietri|1998] l. 

In the models which invoke the global conditions in the ac¬ 
cretion flow to be such that the system finds itself in an unstable 
configuration, the large amplitude fluctuations around the fixed 
point solution are induced. The variability of the disk that re¬ 
flects its global evolution governed by the nonlinear differential 
equations of hydrodynamics, will not be purely stochastic, but 
rather the observed behaviour of the disk will be characterized 
by the deterministic chaos kind of process.This class of models 
which explains the quasi-periodic flares with the global evolu¬ 
tion of the accretion disk, rather than focusing on the flow oscil¬ 
lations in some particular radii, is non-resonant at all. 

The flare-like events of the microquasars such as those ob¬ 
served inIGR 17091-3624 or GRS 1915-1-105 in their’heartbeat’ 
states are very strong, almost coherent low-frequency QPOs. In 
these sources, they are so orderly and quasi-coherent that it is 
possible to see the oscillations even directly in their light curves. 
Generally, it is not possible to identify the QPOs by simply look¬ 
ing at their light curves. A usual method is to apply the Fourier 
transform to the light curve and to study the power spectrum in 
order to identify the QPO as a small peak, or a feature, just up 
the Poisson noise of the power spectrum. Here we apply a novel 
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mathematical method to study the quasi-periodicity and to deter¬ 
mine the nature of the variability of the X-ray sources, i.e., if it 
is stochastic or deterministic. In the latter case, we propose, that 
the underlying mechanism of the deterministic chaotic behaviour 
of the accretion flow is an underlying system of non-linear equa¬ 
tions that governs its global evolution. 


1.1. Accretion disk instabiiity due to the non-iinear 
hydrodynamics 


The problem of chaotic variations in XRBs accretion disks is 
interesting and still unsolved from both theoretical and obser¬ 
vational point of view. First, the classic theory of accretion as 
proposed by Shakura & Sunyaev (|1973|l and|Lightman & Eard 


|ley| \\91A) predicts that the disk is unstable and undergoes the 
limit cycle oscillations, if the viscous stress tensor scales with 
the total (gas plus radiation) pressure, and the global accretion 
rate is large enough for the radiation pressure to dominate. 

In general, the accretion disks may undergo the limit-cycle 
oscillations around a fixed point due to the two main types of 
thermal-viscous instabilities. These are induced either by the 
domination of radiation pressure in the innermost parts of the 
disk, which occurs for high accretion rates, or by the partial ion¬ 
ization of hydrogen in the outer disk parts, for appropriate tem¬ 
peratures. Both these instability types are known for over 40 
years in theoretical astrophysics. 

The hydrogen ionization instability operates in the outer re¬ 
gions of the accretion disk, where the temperatures are in the 
range of log T = 3.5-4 [K] and the hydrogen is partially ion¬ 
ized. Under such conditions the opacities in the plasma depend 
inversely on density and temperature, and non-linear processes 
in the accretion flow are induced. This type of instability leads 
to the outbursts of Dwarf Novae, and Soft X-ray transients, and 
generally its characteristic timescales are on the order of months 
in case of stellar-mass accreting objects (for review see |Lasota| 

( l200T] l). 


Another type of non-linear process, that operates on much 
shorter timescales, is the instability induced by the dominant 
radiation pressure. In classical theory of [Shakura & Sunyaevj 
( 1973| l, the accretion flow structure is based on a prescription 
for the viscous energy dissipation. It assumes that the non-zero 
component of the stress tensor is proportional to the total 
pressure. The latter includes the radiation pressure component 
which scales with temperature as T* and blows up in hot disks 
for large accretion rates. This in turn affects the heating and cool¬ 
ing balance, between the energy dissipation and radiative losses. 
If the accretion rate is small, then most of the disk is gas pressure 
dominated and stable. For large enough accretion rates, there ap¬ 
pears a zone where some of its annuli are dominated by radiation 
pressure and unstable. The larger the global accretion rate, the 
more annuli are affected by the instability, starting at the inner 
edge which is the hottest. 

If there was no stabilizing mechanism, the radiation pressure 
dominated disk would not survive. This is because in such a so¬ 
lution, the decreasing density leads to the runaway temperature 
growth. In consequence, the local accretion rate increases and 
more material is transported inwards. The disk annulus empties 
because of both increasing accretion rate and decreasing density, 
so there is no self-regulation of the disk structure. 

However, the so called ’slim-disk’ solution, where advection 
of energy provides additional source of cooling acts as a stabi¬ 
lizing branch. In this way, advection allows the disk to survive 
and oscillate between the hot and cold states. Such oscillating 
behaviour leads to periodic changes of the disk luminosity. 
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The theoretical calculations of the global, long-term disk 
evolution that are based on the vertically averaged structure, 
clearly show the oscillatory type of behaviour, with character¬ 
istic limit-cycles Paniuk et al.|2002] l. The numerical computa¬ 
tions performed in the 3-D shearing-box simulations with var- 
codes either deny (jKrolik et al. [2007 1 or confirm (Jiang 


[erar][20T3| [Hirose et al.||2009| l that the MHD stress tensor 
scales with the total pressure. Therefore the issue of a-viscosity 
parametrization and contribution of the radiation pressure to the 
energy dissipation is still an open problem. Moreover, the 3-D 
simulations do not tackle the global structure of the flow, and 
therefore are not able to be verified directly with the observa¬ 
tions of astrophysical sources. 


1.2. Observationai support for the accretion instabiiity 


For many years, the only source which undoubtfully showed os¬ 
cillations that are characteristic for the accretion disk instabil¬ 
ity due to the dominant radiation pressure, was the microquasar 
GRS 1915H-105, in some of its observed states ( jBelloni et al.| 
|2000[ jJaniuk et al.||2002| ). The source was thought therefore 
to be unique of its kind, until the very recent discovery of an¬ 
other microquasar, IGR J17091-3624, which was found to be an 
analogue of the previously known microquasar ( [Capitanio et al.| 
2012[ Altamirano & Belloni 2012|l. The recent hydrodynami- 


cal simulations of the global accretion disk evolution confirmed 
that the quasi-periodic flare-like events observed in IGR J17091, 
in some of its variability classes, are also in a good quantitative 


agreement with the radiation pressure instability model (Janiuk 
jet al.|2015» . 


On the other hand, for other sources than these two micro¬ 
quasars, no such detailed analysis was performed, neither obser- 
vationally nor theoretically. However, as proposed by Paniuk] 
|& Czerny|201 l| l, at least eight of the known BH X-ray binaries 
should have their Eddington accretion rates large enough for the 
radiation pressure instability to develop. These other black hole 
binaries can also be promising objects for the radiation pressure 
instability, because they have their Eddington ratio above ~ 0.15. 
The particular value of this critical accretion rate depends on 
whether the viscous torque parametrization is adopted in a mod¬ 
ified version, such as a -sjPg^sPtot, instead of aPtot- Nevertheless, 
for the stellar mass black hole sources, the detailed analysis and 
verification of available observations of various sources, should 
give constraints for the presence of non-linear behaviour of their 
accretion flows. 


Einally, because the instability timescale is scaling with the 
black hole mass, the intermittency in quasars on the timescales of 
hundreds to thousands of years is likely to be of a similar origin 
as in the microquasars. Here also some indirect arguments are 
proposed for the intermittent type of activity in the supermassive 
black holes, e.g. from the studies of the radio maps (Kunert- 
[Bajraszewska et al.|2010[ ). Also, the statistical studies occurred 
to be useful to provide an evidence for the episodic activity in 
AGN ( |C^y et al.|2009| l. 


1.3. Our current approach 

In the current work, we aim to tackle the problem of stochas¬ 
tic versus deterministic nature of the black hole accretion disk 
variability from the analytic and observational point of view. We 
perform the recurrence analysis, which is a powerful tool used to 
study the time series, and known to work in broad range of ap¬ 
plications, ranging from economy to geophysics and medicine 
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( |Eckmann et al.|[l987[ [Marwan et al.||2007) l. The recurrence 
analysis method for processing the time series was used recently 


also for studying chaos in motion of test particles (Kopacek et al. 
|2010[[Sukova|20lT]|Sukova & Semerak|2012[ ). Here, for the first 
time we show that this method has a great power in the studies 
of time series observed from astrophysical sources. 

The essence of recurrence analysis method is to reveal the 
dynamical parameters of the system from the observed time se¬ 
ries. These parameters, such as the Renyi’s entropy (see equa¬ 
tion |A.4| and the following text) or the maximal Lyapunov ex¬ 
ponent, which characterizes the rate of stretching or contracting 
of the attractor (|Ott 


2002 1 , may give indication that the under¬ 


lying variability is deterministic in its nature. This is possible 
if the structure evolution equations of the dynamical system do 
not contain time explicitly, but are nonlinear and have unstable 
steady state solutions. The so-called ’recurrence plot’ contains 
the information about time correlation and has a form of an array 
of points in an NxN square for a time series m,, where i = 

The time series is used to construct the orbits x(i) in the system’s 
c/-dimensional phase space (where x may relate to the systems 
physical state, such as its density or temperature, and d relates 
to the number of non-linear equations that govern its structure 
evolution). The point in the recurrence plot is marked, when¬ 
ever the trajectory returns close to itself, so that x{j) is close to 
jc(0, i.e. closer than some assumed radius e of an t/-dimensional 
sphere. The plot is by definition symmetric with respect to the 
diagonal i = j, and the diagonal lines are not infinite, if only the 
variability is not completely periodic. 

We study the occurrence of long diagonal lines in the recur¬ 
rence plots of observed data series and compare it to the series of 
surrogate data. The latter are produced numerically to have the 
same mean and variance (shuffled surrogates, see below Sec. 
HD or even the same power spectrum (lAAFT surrogates), al¬ 
beit being created by random processes as the permutation of the 
original time series. The surrogate data were used as a counter¬ 
check with the correlation dimension technique and the single¬ 
value decomposition technique already by |Misra et al.| ( |2006| in 
the context of the nonlinear variability in the microquasar GRS 
1915-1-105. Here also, our estimated dynamical invariants, e.g. 
the second order Renyi’s entropy, which is a measure of the pos¬ 
itive Lyapunov exponent and indication of deterministic chaos, 
are determined for both the observed series and the surrogates. 
The example of recurrence plots of several observations and their 


surrogate is given in Fig. A. 3 


We analyze the sample of BHXBs in order to determine the 
nature of their X-ray lightcurves. Our goal is to find the hints 
for deterministic chaotic behaviour of the accreting black hole 
systems, as observed by the advanced X-ray satellites with good 
temporal resolution. Thus we verify the results on the nonlin¬ 
ear behaviour of the GRS 1915-1-105 based on the RXTE data 
( Misra et al.|2006| l and we pursue our analysis to other sources. 
In particular, we study IGR J17091-3624, in the four of its vari¬ 
ability states. We further study the sample of several less certain 
candidates for the globally unstable accretion disk evolution, and 
ultimately we try to answer whether the two microquasars men¬ 
tioned above are unique in their nature. The results of our anal¬ 
ysis should allow for deeper insight into the accretion problem 
and to distinguish whether the variability is governed mainly by 
the set of nonlinear differential equations describing the global 
disk structure, and hence the fundamental physics, or rather by 
the local changes in the flow, and hence environmental condi¬ 
tions different for each source. 

The article is organized as follows. In Section]^ we present 
our full observed sample of black hole X-ray binary sources and 


ObsID 

date 

class 

number 
of points 

01-00 

2.3.2011 

HIMS 

3380 

01-02 

3.3.2011 

SIMS 

5536 

02-00 

4.3.2011 

SIMS/ 

HIMS 

3174 

02-02 

8.3.2011 

HIMS 

6383 

02-03 

10.3.2011 

HIMS 

3305 

03-00 

12.3.2011 

SIMS 

3030 

03-01 

14.3.2011 

IVS 

6682 

04-00 

19.3.2011 

P 

5472 

04-01 

23.3.2011 

IVS 

2297 

04-02 

22.3.2011 

P 

5939 

04-03 

24.3.2011 

IVS 

5499 

05-00 

27.3.2011 

P 

5444 

05-01 

30.3.2011 

P 

2872 

05-02 

25.3.2011 

SIMS 

6620 

05-03 

29.3.2011 

P 

4683 

05-04 

31.3.2011 

P 

5684 

06-00 

2.4.2011 

P 

5227 

06-01 

3.4.2011 

P 

6408 

06-02 

5.4.2011 

P 

5611 

06-03 

6.4.2011 

P 

3617 

07-00 

10.4.2011 

P 

3201 

07-01 

11.4.2011 

P 

6389 

07-02 

12.4.2011 

P 

6388 

08-00 

15.4.2011 

P 

6365 


Table 1. Observations of IGR J17091-3624. The prefix of the ObsID 
is 96420-01-. In the third column the spectral class according to |Pahari| 
|et al.| ( |20T4l l is provided. The number of points of the observation used 
for the analysis in in the last column, the time bin is 0.5s. The data were 
extracted from the generic event mode data using events from PCU2, 
channels 5-24. 


we briefly describe the RXTE data extraction procedure. In Sec- 
tion|^we shortly describe the recurrence analysis method and its 
general properties, while the details about our method are given 
in the Appendix We summarize the main results for the mi¬ 
croquasar IGR J17091-3624 in Sectionj^and in App.[B]we show 
in details how the method was applied to this source. In App. 
we give some more details about the significance of the method, 
depending on the strength of the noise which is always present 
in the observed data sets. In Section]^ we apply the recurrence 
analysis method further to the rest of our observed sample of the 
black hole X-ray binaries. Finally, in Section]^ we discuss our 
results and we present conclusions. 


2. Observations 

We analysed several tens of observations of six black hole X-ray 
binaries, in which the radiation pressure instability is considered 
to develop ( Janiuk & Czerny|2011[ Janiuk et al.|2015| l. We ex¬ 
pect these sources lightcurves to exhibit signatures of non-linear 
behaviour. Below we describe briefly the sources that belong to 
our sample and the method which we used for the extraction of 
the RXTE archival data. 


2.1. Our sample 

The studied sources are the microquasars IGR J17091-3624, 
GRS 1915-H105, GRO 1655-40, XTE J1650-500, XTE J1500- 
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Fig. 1. Lightcurve of the microquasar IGR J17091-3624, for the 
observation ID 96420-01-06-02, extracted in the energy range 2-10 keV 
with time bin dr = 0.5s. We computed the significance of the non-linear 
dynamics detected in this observation .S = 3.72. 
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564 and GX 339-4. Apart from the first one, which went into an 
outburst later in 2011, they were selected from the list given in 
( Janiuk & Czerny |2011[ ) based on the presumably high accretion 
rate to the Eddington rate ratio and high count rate at the same 
time. We included also the source GX 339-4, because higher 
accretion rate was reported by Zdziarski et al.| ( 2004} combined 
with the results of Hynes et al.| ( |2004) l than which was assumed 
in |Janiuk & Czerny ( 201 l| l arid because this source shows also 
other interesting features such as low frequency QPOs with 
changing frequency. 

The first two systems are known to exhibit flare-like events 
in some states (see the references afterwards) and their simi¬ 
lar non-linear behaviour has been studied by different methods 
(|Misra et al.||2004| |2006[ [Altamirano et al.||2011[ [Pahari et'aT] 
2013a|b[). We decided to built-up our method on the sample of 


IGR J17091-3624, whose spectral states have been studied by 
[Pahari et al.| ( |2014] ). The light curve of IGR J17091-3624 during 
the heartbeat state has been modelled by jJaniuk etal. ( 2015| l with 
the model of the disc which undergoes the radiation pressure in¬ 
stability induced oscillations. The wind ejected from the accre¬ 
tion disk regulates the amplitude of the oscillations and increase 
its frequency, so that the simulations are in agreement with the 
observations. Moreover, the stronger winds, which are detected 
in some of the non-heartbeat states via the X-ray spectroscopy, 
suppress the oscillations completely. Therefore, the idea of the 
non-linear instability producing the flares in IGR JI7091-3624 
is strongly supported by this modeling. 

Further, we confirm the results of jMisra et al.j ( |2()04| l for 
some of the observations of GRS 1915-1-105 by our method. 

Later we applied our method for several observations of the 
other four sources. We have not carried out an extended study of 
these sources, we rather aimed to And several examples of obser¬ 
vations, which do or do not show the traces of the non-linear be¬ 
haviour. We took into account the observations with high enough 
count rate and we focused on the rising and declining phases of 
the outbursts. At that time the given source goes through the 
spectral state transition, which can be correlated with an unsta¬ 
ble stage of the accretion disc prone to the internal instability. 

Even though our sample of these four objects is limited and 
serves as a brief example rather than a deep survey, we have 
found indications of non-linearity in several cases, apart from 
the two aforementioned microquasars. 
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2.1.1. IGR J17091-3624 


The microquasar IGR J17091-3624 is a moderately bright tran¬ 
sient X-ray binary (peak flux level at ~20 mCrab in the range 20- 
100 keV) discovered in 2003 ( jKuulkers et al.|2003] l. The 2011 
outburst was the brightest one ever observed from IGR J17091, 
and the source flux increased up to 120 mCrab in the range 2-10 
keV (Capitanio et al. 2012| l. The RXTE/PCA data showed quasi- 
periodic flare-like events occurring at a rate between 25 and 30 
mHz ( [Altamirano et al.|2011| ), which resembled the ‘heartbeat’ 
variation observed previously in the BH binary GRS 1915H-105. 


We tested our method on the sample of observations of 
IGR J17091-3624 selected from those analysed in [Pahari et al.j 
( 2014[ l, which according to their classification belong to differ¬ 
ent spectral classes (hard intermediate state (HIMS), soft inter¬ 
mediate state (SIMS), intermediate variable state (IVS) and vari¬ 
able state (also denoted as p state)). All the observations belong 
into the proposal number 96420, target 01 and they were taken 
by RTXE satellite between March and April 2011. The studied 
observations of IGR 117091-3624 are summarized in Table[T]and 
an exemplary lightcurve of the source is shown in Figure 


2.1.2. GRS 1915-1-105 

GRS 1915H-105is a very well known low mass X-ray binary, 
discovered in 1992 ( [Castro-Tirado et ar][1992| l. Its behaviour 
was interpreted with the time-dependent evolution of an accre¬ 
tion disk, which are thermally and viscously unstable ( |Taam[ 
[et al.|1997| [Janiuk et al.|2000| l. [Belloni et al.[ ( [2000[ ) classified 
for GRS 1915H-105 the 14 classes of variability, out of which 
some are so called ’heartbeat states’, and exhibit regular ampli¬ 
tude periodic oscillations, like the one shown in the lightcurve 
on Figure The flux emission during the heartbeat states of 
the source reaches the Eddington luminosity ( [Done et al.[2004[ l. 
The average mass accretion rate in GRS 1915H-105 on the or¬ 
der of 10** - 3 ■ 10'®g s“' (0.031- 0.63 in Eddington units) was 
estimated by[Neilsen et al.[(|201 l[l. 


2.1.3. GX 339-4 

X-ray nova GX 339-4, discovered in year 1972 ( [Lavrov et al.[ 
1997| l, exhibits high level of X-ray activity. 

The distance was estimated as > 7 kpc by [Zdziarski et al.[ 
( [2004| l, who also reported that the luminosity varied during 
its outbursts (15 outburst from 1987 to 2004) between 0.01 - 
0.25LEdd('7/8kpc)^(10Mo/M). Similar values were also con¬ 
firmed by ( [Kolehmainen et ar[[2011[ l. However, Hynes et al.[ 
P004| l proposed the possibility of the location on the far side of 
the Galaxy with distance > 15 kpc, which would imply that 
the source reaches the Eddington luminosity in the peaks. 


2.1.4. GROJ1655-40 


Low mass X-Ray binary GRO 11655-40 was discovered on Jul, 
27, 1994 with BATSE during the outburst and soon after also its 
optical eclipsing counterpart has been found. [Orosz & Bailyn[ 
\\991} found the mass of the black hole M — 7.02 + 0.22Mq and 
estimated the averaged mass transfer rate as 3.4 ■ lO^^Moyr^* = 
2.16- 10‘^gs-'. 


A strong highly ionized accretion disc wind is reported by 
Miller et al. ( 2008| l on Apr, 1, during the 2005 outburst. The 


origin of the wind is likely not driven by the thermal process. 
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Fig. 2. Lightcurve of the microquasar GRS 1915+105, for the obser¬ 
vation ID 20402-01-03-00, extracted in the full energy range 2-60 keV 
with time bin df = 0.5s from Standard 1 data in thep state (classification 
according to|Misra et al.|2004l>. 



2.1.5. XTEJ1550-564 

XTE J1550-564 was first discovered by the ASM on board 
RXTE on 1998 September 7 during outburst. 

During the 1998 outburst the source exhibited a low and in¬ 
termediate frequency QPOs with evolving frequency in the range 
0.08 - 13 Hz ( |Chakrabarti et al.|20091 l and the accretion rate was 
a significant fraction of Eddington luminosity (|Dunn et al.|2010l 
IHeil et al.|20TT] l. 


2.1.6. XTE J1650-500 


XTE J1650-500 is the black hole X-Ray binary with black hole 
mass4-7.3Mo ( Orosz et al.|2004[|Slany & Stuchlfk 2008[ ). The 
object emits flares with non-thermal energy spectra and occur 
when the persistent luminosity is near 3 ■ 10^^ erg s“^(d/4 kpc)^ 
in 2.8-20 keV range ( [Tomsick et al.|200^ . Homan et al. ( |2006| l 
estimated the distance based on the luminosity during the spec¬ 
tral state transition as 2.6 + 0.7 kpc assuming the black hole mass 
M = 4Mo, for which the Eddington luminosity LEdd ~ 1.1 ■ 10^’ 
erg s“ 


2.2. Data analysis 

We analyzed the online data from the RXTE satellite. We ex¬ 
tracted the time series using Heasoft 6.16 high energy astro¬ 
physics software package. Because we want to compare our 
results for IGR J17091 with the classification of lPahari et al.l 
( |MT4l l, we use the filtering mentioned in that paper using data 
from all xenon layers from the PCU counter number 2 and we 
processed the events from channels 5 to 24 (2-10 keV approx¬ 
imately). Afterwards we selected several observations for the 
other 5 sources, which were observed by RXTE between 1996 
and 2010. These observations were made in different observa¬ 
tional modes with different data available. Eor each source we 
selected data in lower energy band or in full energy range and 
normalized it to count/s/PCU. These sources and details about 
the data extraction are summarized in Table in Section to¬ 
gether with the most important results of our theoretical analysis. 

We adjust the proper binning time to minimize the error and 
simultaneously to not lose the information about oscillations at 
the scale of several seconds. We use the minimal binning time 
between 0.032s and 0.5s. Exact value depends on the flux and 


on the time scale of the variability in the data. Therefore, for 
the flux on the order of hundred counts per second we use the 
minimal binning size 0.5, if the average flux per PCU exceeds 
500 cts/s. Eor 1000 cts/s, we can drop the value to 0.2s, 0.1s or 
even 0.032 s. However, in some cases the variability is slower 
and we can increase the time bin above the minimal level. 

While extracting the data, we use the xdf tool, and Alter the 
data with xtefilt, and we make the ’good time interval’. Eor the 
science array format data files (standard-1 mode or generic sin¬ 
gle bit mode) we extract the data using saextrct command, with 
a proper bin size. The standard-1 mode data are collected in 
the full energy spectrum of PCA (2-60 keV), while the single 
bit mode contains one channel band in the lower energy range 
(it usually contains channels 0-35). Eor the science event format 
data files (generic event data mode), we extracted the lightcurves 
with the seextrct command with the energy range corresponding 
approximately to 2-10 keV. The resulting count rate is normal¬ 
ized to number of PCUs, which were used for the data extraction. 
The data mode used for the extraction is indicated in Table [3 in 
Section|5] 

The background subtraction is ignored in our analysis. This 
is because for the bright sources, it can be safely neglected, and 
also due to the fact that in RXTE observations the background is 
not measured, but simulated numerically. We do not have there¬ 
fore a background measurement in the function of time, while 
the essential part of our analysis is to catch up the changes of 
about 0.5 s time-scales. 

We have not done the spectral analysis of the observations. 
The spectra were studied and in some cases classified by other 
authors. Eor the spectral state we provide in tablej^their results 
or our guess based on the published plots (e.g. on the hardness 
ratio plot from|Miller et al. (|2008]l for the 2005 outburst of GRO 
J1655-40). 


3. Analysis of the chaotic processes 

Our aim is to reveal the important information about the BHXBs 
on the basis of their X-ray light curves using the recurrence 
analysis. 

The recurrence analysis is a well-established method for 
studying the properties of the dynamical systems based on the 
behaviour of its phase space trajectories. More specifically, it 
looks on the “recurrences” in the phase space, which are the mo¬ 
ments, when the phase space trajectory returns close to itself, 
so that the points x(f;) and x(fy) are close. The long diagonal 
lines (of length 1) are formed in the recurrence plot representing 
the recurrence matrix (see relation |A.2| i, when the pairs of suc¬ 
cessive points on the trajectory are close (|| x(f,) - x(fy) ||< e, || 
x(f,+i)-x(fj+i) ||< e,..., II x(f,•+,_!)-x(f/+;_i) ||< e). These lines 
correspond to the case, that the later part of the trajectory evolves 
inside a tube of e diameter along the previous piece. 

In case of a periodic orbit, the recurrence plot would consist 
of infinite diagonals spaced by the period T. Stochastic pro¬ 
cesses yield randomly distributed points, which do not tend to 
form any lines or structures. Chaotic process produces lines with 
finite length, which is related to the maximal Lyapunov expo¬ 
nent, because for chaotic systems nearby trajectories eventually 
diverge outside the e-tube. Therefore, one of the basic but im¬ 
portant quantifiers of the recurrence analysis is the length of the 
longest diagonal line Lmax- We quantify it simply as the number 
of points composing the line. 

The recurrence analysis can be supplied with the time delay 
reconstruction of the phase space trajectory ( |Takens|ri981| l, so 
that the measured time series of one dynamical quantity could 
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be used for studying the underlying dynamical system. From 
that point of view the method is suitable for the analysis of the 
astrophysical data obtained e.g. by X-ray satellites, for which 
only the lightcurve is available. Therefore, the motivation for 
us is also to explore the abilities of recurrence analysis in the 
context of variability of X-ray binaries. 

Here, the underlying dynamical system under study is the ac¬ 
creting gas governed by a global physical parameter, e.g. by the 
mass accretion rate. If the accretion is stationary, the produced 
flux has constant mean and the variability is due to stochastic 
fluctuations. But if the global parameters are such that the intrin¬ 
sic instability of an accretion disk develops, the mass accretion 
rate locally varies in time, which leads to deterministic changes 
of the flux, possibly also overlaid with some stochastic fluctua¬ 
tions. 

Detailed description of the main quantities of the recurrence 
analysis used in this paper is given in Appendix We show 
there several examples of the technique. We also encourage the 
interested reader to look further into the extensive study of the 
method given in |Marwan et al.| ( |2007) l and references therein. 
The details of the method, definitions of the terms and quanti¬ 
ties together with discussion about the possible drawbacks and 
difficulties arising from the presence of noise are given in Ap¬ 
pendices and 

3.1. General approach 

The recurrence analysis technique needs careful setting of the 
parameters. For noisy short data series such estimation is always 
complicated. Therefore instead of the direct estimate of these 
parameters, we rather incorporate the method of surrogate data, 
as proposed by |Theiler et al.|p992[ ). 

We first pose the “null hypothesis” about the measured time 
series, e.g. that the data are product of temporally independent 
white noise or linearly autocorrelated Gaussian noise. Then we 
produce the set of surrogate series, which fulfill this hypothesis, 
but contain as much stochasticity as is possible (at given level 
of the hypothesis), hence there is no further hidden non-linear 
dynamics. 

In the case of the temporally independent white noise, the 
surrogates are created as the random permutation of the values 
of the time series, thus such surrogates (we call them “shuffled 
surrogates”) have the same mean and variance as the original 
time series. In the case of the linearly autocorrelated Gaussian 
noise, the surrogate series are constructed in such a manner, that 
they have the same distribution of values (flux per second per 
PCU for the observed lightcurves) and almost the same spectrum 
( |Theiler et al. |199^ . This can be achieved by an iterative al¬ 
gorithm called Iterative Amplitude Adjusted Fourier Transform 
Algorithm (lAAFT) or its stochastic version SIAAFT, which is 
described e.g. by [Schreiber & Schmitz] ( |2000| l; |Venema et al.| 
( [2006| l. The algorithm iteratively replaces the amplitudes and 
magnitudes of the Fourier coefficients of the surrogate series 
(starting from the white noise) by corresponding values obtained 
from the original time series, so that at the end the surrogate se¬ 
ries is a permutation of the original time series with almost the 
same spectrum. Throughout this paper, we will refer to data se¬ 
ries created by this procedure as simply “surrogates”, or lAAFT 
surrogates. 

For one experimental data series we typically construct one 
hundred of its surrogates. Afterwards, we apply the recurrence 
analysis both on the real experimental lightcurves and on the cor¬ 
responding surrogates and we compare the obtained properties 
for the real and artificial data. 
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If the quantity measured for the real data differs significantly 
from the value obtained for its surrogate, we can reject our initial 
null hypothesis. It means that the non-linear behaviour and per¬ 
haps chaotic nature of the system appears in the case of IA AFT 
surrogates, or non-stochastic but linear behaviour in the case of 
shuffled surrogates. Luckily, there is actually no need for the dis¬ 
criminating criterion to be a particular physical quantity, hence 
the majority of problems with accurate choice of the parameters 
for the analysis are overcome by our approach. The discriminat¬ 
ing criteria are provided by the recurrence analysis of the time 
series. 


3.2. Numerical method 


For some tasks, including the creation of the lAAFT surro¬ 
gates, we employ several proc edures from the publicly available 


software package TISEAN p] (|Hegger et al. 


1999 


Schreiber & 


Schmitz 2000|l. The shuffled surrogates are produced by our 


code written in IDL. The recurrence analysis is performed on the 
basis of the software package described by |Marwan et al.| ( f2007| l; 
Marwan ( |2013 1, which yields also the cumulative histogram of 
diagonal lines. The linear regression for K 2 estimation and other 
post-procession is done by our IDL codes. 

Below, we summarize the main steps of the handling with 
the data. 


- For every observation we extract the lightcurve (discussion 
about the time and energy bin is given in Sect. |2.2| and in 
Table|^. 

- We rescale the extracted light curve to have zero mean and 
unit variance before producing the surrogates for the ease of 
comparison between different observations using procedure 
rescale from TISEAN. 

- We use the procedure mutual from TISEAN in order to find 
appropriate guess of the time delay At = kdt (see App.[A|and 

Fig.IXT}. 

- We use the procedure false_nearest from TISEAN in or¬ 
der to find appropriate guess of the embedding dimension m 
(see App.[A|and Fig. |A.2| i. 

- We produce 100 surrogates (shuffled by our own code, 
lAAFT by the procedure surrogates from TISEAN - or 
both) for each observation. 

- We find the recurrence threshold e^in and e^ax such, that the 
recurrence rate (RR) of the observation is approx. 1% and 
25%, respectivelj0 We create a sequence of e in the range 
(Cmiinfmax) with Constant difference. The number of used 
thresholds ranges between 10 to 40 depending on the length 
of the data series (and thus on CPU time demand of the anal¬ 
ysis). 

- We construct the file with RQA quantifiers and the cumu¬ 
lative histograms of diagonal lines, by the program rp de¬ 
scribed by [Marwan et al.| ( |2007| l; [Marwanj ( |2013| ) with the 
found parameters k, m for all e’s for the observation and the 
surrogates. 

- For each e we compute the estimate of the Renyi’s entropy 

for the observation and for the set of surro¬ 

gates. The Renyi’s entropy of the order a is in general de¬ 
fined as follows: 




1 

I - a 


i=l,N 


( 1 ) 


* http://www.mpipks-dresden.mpg.de/~ tisean/ 

^ The recurrence rate is the ratio of the number of recurrence points to 
all points of the recurrence matrix. 
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where pi is the probability that the random variable has a 
value of i, and basically describes the randomness of the sys¬ 
tem. For detailed explanation of the meaning of entropy K 2 , 
see relations |A.4| and p\.5| in Appendix [A| 

- We compute the average significance .5^ with respect to the 
shuffled surrogates and the averaged significance S with re¬ 
spect to the lAAFT surrogates. The definition of significance 
of chaotic process is taken as; 

- sign(e°^^(e) - (2) 


where Nsi is the number of surrogates, which have only short 
diagonal lines, and is the total number of surrogates, 
gobs gsurr ^j.g jjjg natural logarithms of K 2 entropy for 
the observed and surrogate data, respectively, >Ssi = 3 and 
is the significance computed only from the surrogates, 
which have enough long lines according to the relation 


>Sk2(^) - 


Igobs(^) _ gsu^(g)| 
cr 


( 3 ) 


The significance defined in this way expresses how much the 
value differs from the mean value measured in 

the units of the standard deviation of the set in the 

logarithmic scale crQsutr(f). For further details, the reader is 
asked to read the explanation in App. and the discussion 
before equation[B.2| 


We distinguish between fully stochastic behaviour (those 
observations, which show no significant difference compared 
to the shuffled surrogates, >Sshf < 1.5 or those, which do not 
have enough long lines in the RP), non-stochastic, but linear 
behaviour (observations with a significant result with shuffled 
surrogates, but non-significant result with lAAFT surrogates, 
iSshf > 1.5, < 1.5) and non-linear, possibly chaotic behaviour 

(observations with significant result with respect to the lAAFT 
surrogates, S > 1.5). 

In Tables[^and[^we provide the number of points of the ob¬ 
served time series, which is used for the analysis, the length of 
the longest line for several recurrence rates (the longer the diag¬ 
onal lines are compared to the time delay k, the more regularly 
the system behaves), the averaged significance >Sshf and/or S and 
the number of thresholds for which the average is computed. 


4. Case study of IGR J17091-3624 

We used the sample of observation of IGR J17091-3624 as de¬ 
scribed in Section p.1.1 and in Table[T]for testing the capabilities 
of recurrence analysis to reveal important information about the 
dynamical properties of the source.Details of the procedure are 
given in App. ^and|^ main results are described here and sum¬ 
marized in Table 12] 

We included in our study several observations from the hard 
intermediate state (HIMS), soft intermediate state (SIMS), inter¬ 
mediate variable state (IVS) and variable state (also denoted as 
p state) of the source (see Pahari et ar|p014| l). All observations 
classified by Pahari et al. 1 2014| l as SIMS and HIMS and the ob¬ 
servation 03-01 classified as IVS have low and almost constant 
values of mutual information (see Fig.jATji. Three observations 
(04-01, 04-03 (IVS) and 06-03 (p)) show smooth decrease but 
no second maximum. All other p observations show oscillations 
with period ranging from 10 s to 20 s. 

We found that observations with no oscillations of mutual 
information show no lines long enough for computation of K 2 


ObsID 

N 

5% 

^max 

10% 

^max 

20% 

^max 

S 

Ne 

01-00 

3380 

6 

7 

10 

- 

- 

01-02 

5536 

6 

7 

11 

- 

- 

02-00 

3174 

4 

6 

10 

- 

- 

02-02 

6383 

6 

7 

10 

- 

- 

02-03 

3305 

4 

6 

8 

- 

- 

03-00 

3030 

4 

6 

8 

- 

- 

03-01 

6682 

6 

7 

12 

- 

- 

04-00 

5472 

73 

144 

198 

1.93 

32 

04-01 

2297 

6 

8 

13 

- 

- 

04-02 

5939 

33 

117 

274 

3.21 

27 

04-03 

5499 

6 

10 

14 

- 

- 

05-00 

5444 

12 

22 

100 

2.35 

18 

05-01 

2872 

10 

26 

89 

1.08 

15 

05-02 

6620 

6 

6 

10 

- 

- 

05-03 

4683 

13 

72 

126 

2.45 

19 

05-04 

5684 

13 

68 

164 

3.18 

19 

06-00 

5227 

43 

130 

215 

2.73 

28 

06-01 

6408 

125 

280 

344 

3.29 

35 

06-02 

5611 

134 

207 

391 

3.72 

35 

06-03 

3617 

19 

77 

79 

1.21 

18 

07-00 

3201 

137 

223 

290 

2.26 

34 

07-01 

6389 

125 

171 

289 

3.04 

34 

07-02 

6388 

106 

180 

285 

3.13 

36 

08-00 

6365 

35 

51 

100 

1.63 

17 


Table 2. Results of our analysis for IGR J17091-3624. The prefix of 
the ObsID is 96420-01-. We used the RP parameters m = 10, At = 7s 
for every observation. The number of points N of the observation used 
for the analysis is in the second column, the time bin is 0.5s. Next, we 
present the number of points composing the longest diagonal line for 
three different thresholds yielding the recurrence rate approx. 5% (third 
column), 10% (fourth column) and 20% (fifth column). The average 
significance S and the number of different thresholds M used for the 
averaging is in the two last columns. 


(the criteria for the length and number of lines are described 
in App. 1^ and their L^ax behaves similarly like that of both 
types of surrogates.. Therefore, these observations are consis¬ 
tent with the hypothesis of temporally independent white noise 
(see Fig. |A.6| l and the source did not show any dynamical evolu¬ 
tion during that time. 


The observations 06-03 and 05-01 do contain enough long 
lines, so that the significance can be computed. In both cases 
non-significant results are obtained with respect to the lAAFT 
surrogates. However, we obtainsignificance with respect to shuf¬ 
fled surrogates >Sshf = 3 for both,computing the average using 
Vf = 13 and = 15 different thresholds for 06-03 and 05-01, 
respectively. Thus, we classify these two observations as a non¬ 
stochastic linear process. As mentioned in Appendix ^ this can 
also mean, that these observations refer to the period, when the 
source was following regular or “sticky trajectory” of the gener¬ 
ally non-linear dynamical system. Such orbit for very long time 
evolves very close to a regular island, sharing its main dynamical 
properties, but eventually departures to another part of the phase 
space, differing significantly from the regular motion (for more 
details see e.g. Semerak & Sukova (2012 1 ). Pahari et al. ( |2014| l 
classified these observations as the variable p state. 


The rest of the p state observations show significant result 
with S > 1.5. The values range between 1.63 for ObsID 08-00 
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Fig. 3. Lightcurve of the source GX 339-4. The observation ID is 
95409-01-16-05 and the source was in its soft-intermediate state. Data 
were extracted in the energy range 2-10 keV (channels 5-25). We com¬ 
puted the significance of the non-linear dynamics detected in this obser¬ 
vation S = 2.90. 


Fig. 4. Lightcurve of the source XTE 11550-564. The observation 
ID is 30188-06-03-00 and the data were extracted in the full energy 
range from standard-1 mode data. We computed the significance of the 
non-linear dynamics detected in this observation S = 2.89. 


to Ji.ll for ObsID 06-02. For these observations we claim that 
the source was showing non-linear behaviour. 

For illustration we also provide the information about the 
length of the longest diagonal line for different values of RR 
in Table 1^ Higher ratio Lmax/^ means, that the system evolves 
similarly for longer time compared to the period of the oscilla¬ 
tions (i.e. that the phase trajectory in two different times evolves 
within e-tube). Further details are given in App. 


5. Using the method on other sources 

The last step of our study is to apply our method to other sources 
listed in Section |2] We have chosen several observations of five 
objects, GRS 1915-1-105, GRO J1655-40, XTE J1650-500, XTE 
J1550-564 and GX 339-4, which are summarized in Table|^ 

At first, we examined observations of GRS 19154-105, which 
is a well studied source and several papers about the properties 
of the lightcurves with respect to the non-linear nature of the 
source have been published. We picked the observations studied 
by Misra et al. ( 2004| l, who classified them according to their 
correlation dimension into three different groups, namely those 
which show the chaotic behaviour (spectral classes /(, A, k, and 
//), non-stochastic behaviour (spectral classes 9, p, v, a, and 6), 
and stochastic behaviour (spectral classes <&, y and x)- An ex¬ 
emplary lightcurve taken in the class p is shown in Eigure]^ 

In agreement with their results, we found > 1.5 for all ob¬ 
servations belonging into the first and second group and we have 
not obtained diagonal lines long enough for the K 2 estimation 
for observation from class <&. The averaged significance ranges 
between 1.6 and 5.1 for different observations. The variability 
of GRS 1915-H105 is slower than in case of IGR J17091-3624, 
so that the time delay At - kdfis longer and ranges between 10s 
to even 125s. Therefore, the ratio of the number of points of 
the data set and the embedding delay N/k is lower than for IGR 
117091-3624 for similar length of observation (lower number of 
“cycles” is seen), hence lower significance could be expected 
for the same type of behaviour. The relatively high values of av¬ 
eraged significance thus provides quite strong evidence for the 
non-linear behaviour of the source. 

Eor the black hole candidate GX 339-4 we analysed four ob¬ 
servations from four different spectral classes classified accord¬ 


ing to Nandi et al. ( 2012| l. In Eigurej^we show the exemplary 
lightcurve, as taken on Sep. 24, 2010, when the source was in 
its soft-intermediate state. Here we have found the averaged sig¬ 
nificance equal to 2.90. Eor the observation taken on Mar. 26, 
2010, the averaged significance is equal to 2.19. The third obser¬ 
vation yields = 1.59, which is close to our chosen threshold 
for the non-linear behaviour. We have also computed the signif¬ 
icance with respect to the shuffled surrogates .Sshf, which show 
non-stochastic behaviour in all three cases. Eor the last observa¬ 
tion of GX 339-4, non-significant result is obtained with respect 
to both types of surrogates, hence the variability of the source in 
this state is of stochastic origin. 

Figure 1^ shows the lightcurve of the source XTE J1550- 
564, as taken on Sep 08, 1998. The source has been previ¬ 
ously studied by [Sobczak et al. ( 2000[ ), and the lightcurve of 
observation 30191-01-14-00 shown therein looks similar to the 
variable state of IGR J17091-3624, however on much smaller 
time scales. Thus, we needed to extract this lightcurve with very 
small time bin (df = 0.032s), which increases the influence of 
noise. On the other hand, this data set has much higher num¬ 
ber of points. Our analysis of this observation (with the highest 
N/k ratio) yields the averaged significance S = 6.14, which is 
the highestamong the sample. Eor comparison, we also provide 
the significance with respect to the shuffled surrogates, which is 
iSshf = 25.17. Also three other observations yield quite high 
significance around 2.7. Eor the lightcurve presented in Eig- 
ure the computed significance of the non-linear dynamics is 
S = 2.89. 

The other two observations of XTE J1550-564 show non¬ 
significant results (S - 1.00 and S = 0.30). However, their RPs 
contain quite long lines compared to their shuffled surrogates, 
which do not show any lines longer than At with the same recur¬ 
rence parameters {Nsf; = 0 for all e), so that >Sshf = 3.00 in both 
cases. This indicates the non-stochastic behaviour of the source 
during these states. 

Eigure shows an exemplary lightcurve of the source GRO 
11655-40, as taken by RXTE on Aug 01, 1996 (10255-01-04- 
OOB). For this data set, the significance of the non-linear dynam¬ 
ical process that governs its variability, is above S = 4.9. How¬ 
ever, in this case the high value of significance could partially be 
caused by the fact, that during this observation several episodes 
with different behaviour occur, so that the evolution is not sta- 
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source 

ObsID 

date 

data 

mode 

state/ 

class 

N 

df [s] 

m 

k 

20% 

^max 

S 

*^shf 

N, 

GRS 

10408-01-08-00 

21.05.1996 

Stdl 


289 

5 

8 

12 

438 

3.16 

- 

15 

GRS 

10408-01-10-00B 

26.05.1996 

SB,ch0-35 


6585 

0.5 

6 

100 

1439 

5.00 

- 

33 

GRS 

10408-01-10-00C 

26.05.1996 

SB,ch0-35 

P 

6547 

0.5 

6 

100 

1685 

2.69 

- 

33 

GRS 

10408-01-12-00B 

05.06.1996 

Stdl 

O 

5466 

0.5 

8 

50 

33 

- 

- 

- 

GRS 

10408-01-12-00C 

05.06.1996 

Stdl 

O 

6485 

0.5 

8 

50 

42 

- 

- 

- 

GRS 

10408-01-12-00D 

05.06.1996 

Stdl 

o 

6626 

0.5 

8 

50 

33 

- 

- 

- 

GRS 

10408-01-17-00 

22.06.1996 

Stdl 

6 

6812 

0.5 

8 

60 

300 

3.17 

- 

17 

GRS 

20402-01-03-00B 

19.11.1996 

Stdl 

P 

6543 

0.5 

8 

18 

458 

1.63 

- 

33 

GRS 

20402-01-03-00C 

19.11.1996 

Stdl 

P 

3720 

0.5 

8 

18 

372 

2.47 

- 

34 

GRS 

20402-01-03-00E 

19.11.1996 

Stdl 

P 

3407 

0.5 

8 

18 

309 

2.68 

- 

30 

GRS 

20402-01-33-00A 

18.06.1997 

SB,ch0-35 

K 

6238 

0.5 

6 

50 

480 

1.63 

- 

27 

GRS 

20402-01-33-00B 

18.06.1997 

SB,ch0-35 

K 

6690 

0.5 

6 

50 

554 

1.60 

- 

27 

GRS 

20402-01-37-01A 

12.07.1997 

SB,ch0-35 

A 

6619 

0.5 

6 

250 

1474 

2.56 

- 

14 

GRS 

20402-01-37-0 IB 

12.07.1997 

SB,ch0-35 

A 

6615 

0.5 

6 

250 

900 

2.80 

- 

15 

GRS 

20402-01-45-02B 

05.09.1997 

SB,ch0-35 

0 

6612 

0.5 

6 

100 

985 

3.29 

- 

27 

GRS 

20402-01-45-02C 

05.09.1997 

SB,ch0-35 

0 

6612 

0.5 

6 

100 

916 

5.09 

- 

30 

GRS 

20402-01-45-02D 

05.09.1997 

SB,ch0-35 

0 

6466 

0.5 

6 

100 

795 

4.10 

- 

23 

GX 

95409-01-12-00 

26.03.2010 

GE,ch5-25 

HS 

3170 

0.5 

8 

5 

89 

2.19 

5.38 

22 

GX 

95409-01-14-06 

13.04.2010 

GE,ch5-25 

HIMS 

3155 

0.5 

8 

5 

26 

0.52 

0.75 

8 

GX 

95409-01-16-05 

29.04.2010 

GE,ch5-25 

SIMS 

1625 

1.5 

8 

12 

66 

2.90 

3.00 

14 

GX 

95409-01-35-02 

19.09.2010 

GE,ch5-25 

SS 

2906 

0.5 

8 

5 

42 

1.59 

2.05 

8 

XTE 1 

30188-06-02-00A 

07.09.1998 

Stdl 

HS* 

3589 

0.125 

7 

15 

92 

1.00 

3.00 

8 

XTE 1 

30188-06-02-00B 

07.09.1998 

Stdl 

HS* 

5019 

0.125 

7 

15 

107 

0.30 

3.00 

15 

XTE 1 

30188-06-03-00A 

08.09.1998 

Stdl 

HS* 

7102 

0.25 

8 

7 

115 

2.89 

- 

25 

XTE 1 

30188-06-03-00C 

08.09.1998 

Stdl 

HS* 

13595 

0.25 

8 

7 

125 

2.76 

- 

25 

XTE 1 

30188-06-01-01 

09.09.1998 

Stdl 

HS* 

10718 

0.125 

8 

4 

114 

2.52 

- 

22 

XTE 1 

30191-01-14-00 

28.09.1998 

SB,chO-17 

VHS 

109388 

0.032 

9 

3 

151 

6.14 

25.17 

13 

XTE 2 

60113-01-13-01 

19.09.2001 

Stdl 

HS 

5463 

0.125 

8 

4 

38 

0.14 

0.29 

15 

XTE 2 

60113-01-13-02 

19.09.2001 

Stdl 

HS 

7127 

0.125 

8 

4 

42 

0.17 

0.33 

18 

XTE 2 

60113-01-18-00 

24.09.2001 

GE,ch5-24 

H/ST 

17715 

0.1 

8 

4 

19 

-0.30 

0.22 

10 

XTE 2 

60113-01-18-01A 

24.09.2001 

GE,ch5-24 

H/ST 

19803 

0.1 

8 

4 

52 

0.72 

1.20 

10 

XTE 2 

60113-01-18-01B 

24.09.2001 

GE,ch5-24 

H/ST 

15591 

0.1 

8 

4 

59 

0.34 

0.74 

9 

XTE 2 

60113-01-18-02 

24.09.2001 

GE,ch5-24 

H/ST 

8581 

0.1 

8 

4 

46 

0.93 

1.14 

10 

XTE 2 

60113-01-19-00 

25.09.2001 

GE,ch5-24 

H/ST 

19191 

0.1 

8 

4 

44 

-0.37 

-0.17 

10 

XTE 2 

60113-19-01 

26.09.2001 

Stdl 

SS 

1713 

0.25 

8 

4 

31 

0.13 

0.39 

12 

XTE 2 

60113-19-04 

27.09.2001 

Stdl 

SS 

1346 

0.25 

8 

4 

22 

-0.28 

-0.25 

11 

XTE 2 

60113-29-00 

17.10.2001 

Stdl 

SS 

6413 

0.25 

8 

5 

41 

0.42 

0.43 

14 

GRO 

10255-01-04-00A 

01.08.1996 

SB,ch0-35 

VHS 

3287 

0.5 

7 

10 

1602 

2.21 

- 

24 

GRO 

10255-01-04-00B 

01.08.1996 

SB,ch0-35 

VHS 

6125 

0.5 

7 

10 

962 

4.97 

- 

24 

GRO 

10255-01-04-00C 

01.08.1996 

SB,ch0-35 

VHS 

6837 

0.5 

7 

10 

201 

-0.45 

3.00 

18 

GRO 

10255-01-18-00A 

02.11.1996 

SB,ch0-35 

VHS 

642 

5 

7 

5 

58 

1.85 

- 

14 

GRO 

10255-01-18-00B 

02.11.1996 

SB,ch0-35 

VHS 

5663 

0.5 

7 

9 

65 

2.05 

- 

11 

GRO 

20402-02-14-OOA 

28.05.1997 

SB,ch0-35 

HSS 

699 

5 

7 

10 

52 

1.35 

3.00 

12 

GRO 

20402-02-14-OOB 

28.05.1997 

SB,ch0-35 

HSS 

620 

5 

7 

10 

49 

-0.10 

3.00 

11 

GRO 

90058-16-02-00 

21.02.2005 

Stdl 

HS* 

6502 

0.5 

7 

10 

17 

0.18 

1.98 

3 

GRO 

90058-16-03-00A 

23.02.2005 

GE,ch5-24 

HS* 

4534 

0.5 

7 

5 

32 

0.22 

1.74 

24 

GRO 

90058-16-03-00B 

23.02.2005 

GE,ch5-24 

HS* 

1973 

0.5 

7 

5 

27 

0.46 

1.57 

17 

GRO 

90058-16-07-00A 

26.02.2005 

GE,ch5-24 

HS* 

4111 

0.5 

7 

5 

36 

0.17 

-0.18 

23 

GRO 

90058-16-07-00B 

26.02.2005 

GE,ch5-24 

HS* 

2388 

0.5 

7 

5 

38 

1.86 

- 

21 

GRO 

90019-02-01-00 

13.03.2005 

Stdl 

SIMS* 

1549 

0.5 

6 

12 

45 

2.33 

- 

5 


Table 3. Table of the RXTE observations of other sources, the meaning of the shortcuts is the following: GRS = GRS 1915+105, GX = GX 339-4, 
XTE 1 = XTE J1550-564, XTE 2 = XTE J1650-500, GRO = GRO J1655-40, Stdl = Standard 1 data mode, SB,ch0-35 = Single bit data mode 
containing energy channels 0-35 and GE,ch5-25 = generic event data mode using energy channels 0-25 from PCU2. The classification of the 
source GX 339-4 is given according to |Nandi et al.|p012^ (HS = hard state, HIMS = hard-intermediate state, SIMS = soft-intermediate state, 
SS = soft state), the classification of the source GRS 1915+105 is given according to |Misra et al.| ( |20d4^ , the classification of GRO J1655-40 
is according to |Sobczak et al.| l |1999a[ l (VHS = very high state with flaring, HSS = high/soft state) for the 1996/1997 measurement. During the 
2005 outburst we used the overall plot of the hardness ratio given in [Miller et al.| ( (2008] l for an estimation of the spectral state (HS* - hard state, 
SIMS* - soft-intermediate state). Similarly we used the classification, plots and discussion given by |Sobczak et al.| ( |1999b^ for XTE J1550-564 
and by |Homan et al.|f2003^ for XTE J1650-500 (H/ST - hard to soft state transition). In the 6-th and 7-th columns the used number of points of 
the lightcurve N with the time bin df are summarized. In the next columns,_the parameters (m,k) od the recurrence analysis and its results - the 
length of the longest diagonal line for RR~20%, the averaged significance <S and the averaged significance with respect to the shuffled surrogates 
.Sshf computed for different values of e - are presented. 
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Fig. 5. Lightcurve of the source GRO 11655-40. The observation ID 
is 10255-01-04-00 and the data were extracted in the energy range 2-13 
keV (channels 5-35). We computed the significance of the non-linear 
dynamics detected in this observation S = 4.97. 


Fig. 6. Lightcurve of the source XTE 11650-500. The observation ID 
is 60113-01-18-02 and the data were extracted in the energy range 2-10 
keV (channels 5-24). We computed the significance of the non-linear 
dynamics detected in this observation S = 0.93. 


tionary. In case of non-stationary data, the theoretical problem 
is much more complicated and using the surrogates created with 
the same spectrum may not be appropriate ( |Theiler et al.|1992| l. 

Altogether, for this source we have found six observations 
with > 1.5 and seven other with non-significant results. For 
the latter group, we computed also the shuffled significance. 
The obsID 90058-16-07-00A shows again non-significant re¬ 
sult, hence it is of stochastic origin. Three other observa¬ 
tions (90058-16-03-00A, 90058-16-03-00B and 90058-16-02- 
00) yield moderate results <Sshf e (1.5,2.0), indicating a weakly 
non-stochastic origin. The shuffled surrogates of the last three 
observations (10255-01-04-00C, 20402-02-14-00A and 20402- 
02-14-00B) have no long lines, leading to .Sshf = 3.0, which is a 
strong evidence of a non-stochastic behaviour of the source. 

Finally, within the studied sample of observations of XTE 
J1650-500, no significant traces of non-linear behaviour have 
been found. Figure shows one of its exemplary lightcurves, 
from the observation ID 60113-01-08-02, taken on Sep 24, 2001, 
for which the computed significance of chaotm process was the 
highest. In any case, it was always below S - 1.0. For this 
source also the significance with respect to the shuffled surro¬ 
gates is always below <Sshf =1.5. Therefore we have not found 
any observation of this source, which departures from stochastic 
behaviour. 


6. Discussion and conciusions 

6.1. Discussion 

In this article, we used the recurrence analysis method to study 
the non-linear behaviour of several X-ray sources. They are the 
black hole candidates, in which their X-ray luminosity origi¬ 
nates in an accretion disk. The latter may be responsible for a 
short-timescale, quasi-regular oscillations of the emitted energy 
flux, through some physical process related to its hydrodynam- 
ical evolution. The possible trigger for non-linear behaviour 
could be the moment, when the disc enters the parameter region, 
where the thermal-viscous instabilities occur. These instabilities 
lead to global limit cycle oscillations of the disc configuration 
(including density, mass accretion rate, temperature, etc.). Such 
unstable evolution of the disc translates into non-linear patterns 
in the light curve. 
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Two microquasars, GRS 19154-105, and recently discovered 
IGR117091-3624, are already firm candidates for the limit-cycle 
oscillations that are called also the ’heartbeat’ state (see e.g., 
Neilsen et al. |201 l|l. Our present analysis confirmed that the 


variability in these sources is significantly governed by the non¬ 
linear dynamics of accretion process. Other sources of that kind, 
which accrete at presumably high accretion rates and exhibit 
soft, disk-dominated, states, have not yet been studied exten¬ 
sively enough to tackle this problem quantitatively. The ques¬ 
tion however remained, whether GRS 1915-1-105, and later its 
analogue, IGR 117091-3624, are the only two sources of this 
kind, or rather such behaviour can be common to all accreting 
black holes, if only the physical conditions which trigger the 
non-linearity, are similar in there. 

For other X-ray binaries than these two microquasars, the sit¬ 
uation is not so obvious. However, |laniuk & Czerny| ( |201 1[ ) dis¬ 
cussed the theoretical aspects of the thermal-viscous instabilities 
in the accretion disks and presented observational constraints for 
a sample of Galactic black hole systems. In particular, for the 
radiation pressure instability, a few other sources were hinted, 
based on their estimated Eddington ratios and observed prop¬ 
erties. It should be stressed, that the accretion rate at exactly 
Eddington rate is not required by the theory for the radiation 
pressure instability to develop. Therefore we aim to search for 
the non-linear hydrodynamics of the accretion disk hidden in the 
data also for these other sources. 

The black hole X-ray transient source XTE J1650-500 has 
shown its only outburst in the year 2001, and there were not 
many observations available for that source. Homan et al.| ( |2003) 
detected the high-frequency variability in XTE J1650-500 and 
interpreted them as the orbital frequency at the innermost stable 
orbit around a Schwarzschild black hole, while [Kalemci et al.| 
( |2003| l suggested that the accretion flow geometry is that of the 
disk-jet-corona system. There is not Arm estimation in the litera¬ 
ture about the value of the accretion rate in this source. ITomsickI 


et al. (20031 detected in this source the X-ray flares which have 


durations between 62 and 215 s and peak fluxes that are 5-24 
times higher than the persistent flux. Therefore |Janiuk & Czerny | 
( |2011| l tentatively hinted at this X-ray binary as a prospective 
candidate for the disk instability induced oscillations, because 
of these timescales. These flares have however non-thermal 
energy spectra, and possibly do not originate in the accretion 
disk/corona at all. They also occur at luminosities about thou- 
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sand times lower than the peak outburst luminosity, probably at 
the level about 2 ■ lO^^LEdd- Therefore they are not good candi¬ 
dates for the radiation pressure instability. Tomsick et al. P003| l 
also reported about aperiodic oscillation with 14 days time scale. 
They related these variations of X-ray light curve with the ion¬ 
ization instability. However such a long variability could not be 
studied by our present method, because no continuous observa¬ 
tion of required length (at least one order of magnitude longer 
than the period) is available. Our current analysis was based on 
the earlier observations, from up to October 2001, with a higher 
count rate. It showed that in this source the significance of the 
underlying non-linear dynamics is very low. We would therefore 
rather reject it as a candidate for the limit-cycle type of oscilla¬ 
tions. This source requires further monitoring to put more firm 
conclusions about the nature of its variability. 

For the other sources in our study, the significance of the 
underlying deterministic chaotic process, being the intrinsic to 
the accretion disk reason for its observed variability, is much 
higher. This result is in full agreement with the observational 
facts, noted during the time of over 20 years of collecting data 
by the X-ray satellites. 

In XTE J1550-564, at the peak of the outburst the luminosity 
is close to LEdd ( |Heil et al.|201 l| l, and the quasi-periodic oscilla¬ 
tion (QPO) observed in this source during its 1998 outburst were 
found to be linked to the disc count rate ( |Sobczak et al.|2000l l. 
XTE J1550-564 is classified as a microquasar on the basis of 
its large-scaled moving jets, detected a t X-ray and radio (|Corbel| 
et al.|200^ i, similarly with GX 339-4 (jCorbel & Eender|2002|l. 


We expect that these two objects should have similar character¬ 
istics of the disk variability, and along with the two well stud¬ 
ied microquasars the non-linear dynamical processes in these 
sources should also occur. Our current analysis does confirm 
these expectations. 

The spectral state of GRO J1655-40 was initially not deter¬ 
mined (i.e., the low, high, or very high state) for its first outburst 


in 1994, based on the BATSE observation (Crary et al. 


19961. 


However, [Remillard et al.| ( | 1999) 1 after the analysis of the data 
spanning the X-ray outburst from 1996-1997, estimated the un¬ 
absorbed X-ray luminosity to above ~ 0.2LEdd during the peak. 
Such luminosity and accretion rate is quite sufficient for the in¬ 
stability of the accretion disk, which develops within the inner¬ 
most 60-80 Schwarzschild radii from the black hole, and leads to 
moderate luminosity oscillations ( jJaniuk & Czerny |201 l) l. |Men-| 
|dez et al.| ( 1998| l studied the evolution of GRO J1655-40 through 
the high, intermediate, and low state. This source at the begin¬ 
ning of its decay might have even shown signatures of a very 
high state, just like other black hole candidates. |Sobczak et al.| 
( |1999a| l presented the full spectral analysis of the RXTE data for 
1996-1997 outburst of GRO J1655-40 and showed that during 
the high/soft state, its spectrum is dominated by the soft ther¬ 
mal emission from the accretion disk. Comparing the two above 
mentioned black hole binaries, [Sobczak et al.| ( |2000) l found that 
in both sources the QPO frequency is correlated with the disk 
flux, and hence with the rate of mass accretion through the disk, 
confirming therefore that the intrinsic mechanism of this vari¬ 
ability is related to the accretion process. Most recently, |Utt-| 
ley & Klein-Wolt|(20151 reported about an unusual soft state of 


GRO J1655-40, observed during its 2005 outburst by the RXTE. 
Chandra X-ray grating observations have revealed a high mass- 
outflow accretion disc wind in this state, which can be related 
to the accretion disk (in)stability evolution Paniuk et al.|2015y 
Indeed, |Miller et al.| ( |2008| l showed, that during its 2005 out¬ 
burst, GRO J1655-40 ejected massive winds, possibly driven by 
magnetic process in the accretion disk. 


The fact, that we have found non-linear dynamics hidden in 
the light curves of almost all the studied sources means, that the 
evolution of the accretion disc, and possibly corona, is important 
for the outgoing radiation and that the accretion flows is influ¬ 
enced by a physical instability, at least at some specific states. 

In App. [C]we summarized the results presented in 
|& Janiuk|2015, in preparation] l, where we test the capabilities of 
our method on two numerical trajectories. We showed how the 
method works with simulated trajectories of a complicated non¬ 
linear system (motion of the test particle in the held of a black 
hole surrounded by thin massive disc given by exact solution of 
Einstein equations), from which one is a regular orbit, while the 
other one is chaotic. The chaotic orbit shows higher significance 
of non-linear dynamics. We studied the influence of noise on the 
output of the method. While the significance for regular orbit 
quickly drops below unity, the significance of chaotic motion 
remains high up to comparable variance of the noise and the data. 


(Sukova 


6.2. Conclusions 


We applied the recurrence analysis to the observations of six 
black hole X-ray binaries detected by RXTE satellite. 

We developed a method for distinguishing between stochas¬ 
tic, non-stochastic linear and non-linear processes using the 
comparison of the quantification of the recurrence matrices 
for the real and surrogate data. 

We tested our method on the sample of observations of the 
microquasa r IGR J17091-3624 . which spectral states were 


provided by Pahari et al. |2014 1 . Results with > 1.5 for the 


“heartbeat” variable p state were obtained, which indicates 
the non-linear dynamics. Hence, our results also corroborate 
the validity of simulations of radiation pressure instability in¬ 
duced oscillation of the accretion disc in this source reported 
by pamukeTaL] (|20T5]| . 

We examined several observations of other five micro¬ 
quasars. Aside from the well-studied binary GRS 1915-1-105, 
we found significant traces of non-linear dynamics also in 
three sources (GX 339-4, XTE J1550-564 and GRO J1655- 
40). 

The non-linear behaviour of the lightcurve during some ob¬ 
servations gives the evidence, that the temporal evolution 
of the accretion flow in the binaries is governed by low¬ 
dimensional system of non-linear equations and the dynam¬ 
ics is chaotic. Hence the variability in the lightcurve is driven 
by a dynamical system of the accreting gas described by 
a few parameters (e.g., the accretion rate) rather than by 
stochastic variations coming from random processes (flares) 
in the disc. Possible explanation is that the accretion disc 
is in the state prone to the thermal-viscous instability and is 
undergoing the in duced limit cycle osc illations. 

As pointed out in Misra et al. (20041 the knowledge about 
the chaotic nature of the system can enlighten the fundamen¬ 
tals of the physical processes yielding the observed radiation. 
The rough estimate of the Renyi’s entropy, which is related 
to the maximal Lyapunov exponent, can be compared with 
the estimates for magnetohydrodynamical simulations of the 
accretion flow to further confirm or refute the accretion mod¬ 
els. 


Our analysis of course does not answer the question of how 
the appearance or absence of non-linear chaotic variability is re¬ 
lated to the evolution of a given X-ray binary during its outburst 
and the spectral state transitions. Such analysis is extremely 
complex and currently beyond the scope of this article, but defi¬ 
nitely worth further studies. Eor now, our tentative picture would 
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be that the purely stochastic variability occurs when the source 
is in its hard state, at the very beginning of the outburst. In the 
disk dominated soft state, as well as in the intermediate states, 
the non-linear variability due to chaotic process in the accretion 
disk may appear. However, even in this state, the intrinsic coro¬ 
nal emission, which should be rather stochastic, as well as pos¬ 
sibly reprocession of the hard X-ray flux in the disk, may sig¬ 
nificantly affect the observed flux and hence our results. On the 
other hand, if the corona above the disk forms through the evapo¬ 
ration of the upper layers of an accretion disk, then its variability 
will be related to the disk variability, with a possible short time 
delay ( |Janiuk & Czerny|20051 l. 
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Appendix A: Recurrence analysis of the observed 
time series of IGR J17091-3624 

The recurrence analysis is based on the study of the times, when 
the trajectory returns close to itself in the phase space (closer 
than a certain recurrence threshold e). However, the trajectory 
has to be reconstructed from the observed time series with the 
time delay technique. Hence the resulting phase space vector is 
given as 

y(t) - {x{t), x{t + At), x(t + 2At),... ,x{t + im - l)Af)), (A.l) 


where x{t) is the time series. At is the embedding delay and m is 
the embedding dimension. When working with measured time 
series, the embedding delay At has to be an integer multiple of 
the binning time of the data df, hence At - k df. The determi¬ 
nation of these parameters is a subtle issue ( Kennel et al.|1992| 
|Cao||1997l [Thiel et ar]|2004[ [Marwan et al.||2007| ). The usual 
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Fig. A.l. The dependence of mutual information on the time delay 
At for the set of observations of IGR J17091-3624 listed in Table [f| 
Observations with low and constant mutual information are plotted with 
black dashed line, observations with oscillating mutual information are 
plotted with grey solid line and those with intermediate behaviour are 
plotted by yellow dotted line. The vertical line indicates the chosen time 
delay for our analysis. 


Fig. A.4. The cumulative histogram pde, 1) for the observation 06-02 
computed for e = 2.8. The length of the lines / is given in number of 
points of the lines (thus the total time of the line is given as t; = I At). 
The estimate of the Renyi’s entropy K 2 = 0.0699 ± 0.0003 is obtained 
as the linear regression of the central part of the cumulative histogram 
(indicated by solid red line, the dashed red lines mark the used part of 
the histogram). 



Fig. A.2. The dependence of the ratio of false nearest neighbours on 
the embedding dimension m for the set of observations of IGR J17091- 
3624 for the choice R,oi = 10. 


practise is to determine the time delay as the first distinct mini¬ 
mum of the time delayed mutual information, for which we use 
the procedure mutual from the TISEAN package. In Fig. |A.1| 
the dependence of the mutual information on the time delay At 
is given for the set of observations of IGR J17091-3624. The 
observations show two kinds of behaviour. Mutual information 
of some lightcurves is very low even for small time delays and is 
roughly constant for all values of At (black dotted lines), whereas 
the course of mutual information of others is smooth and shows 
oscillations with period ranging from 10 s to 20 s (grey solid 
lines). Three observations show intermediate behaviour, there 
is a smooth decrease of the mutual information, but no clear sec¬ 
ond maximum, the values then rather stabilize around a constant 
value of the first minimum (yellow dotted line). Because we 
want to compare the results of the analysis with the same pa¬ 
rameters for all observations, we choose the time delay to be 
At - 7 s (for df = 0.5 s is k = 14) as a reasonable value for each 
lightcurve. 

For finding the appropriate embedding dimension, the well 
known method is to determine the number of false nearest neigh¬ 


bours ( jKennel et al.|1992) l. This procedure is based on the fact, 
that for too low embedding dimension the points, which are not 
close on the attractor, can be projected into the smaller space 
as very close neighbours. When the dimension is increased, the 
ambiguity of their mutual position is removed and their distance 
grows substantially. Hence, we set a tolerance threshold Rtoi 
and denote the point as false neighbour, if the distance in m and 
m+1 dimension increase by a factor bigger than Rtoi ■ In this way 
we investigate the change of distance for each point of the trajec¬ 
tory and its nearest neighbour in the reconstructed phase space. 
If the dimension is big enough to cover the whole attractor, the 
fraction of false neighbours goes to zero. We use the procedure 
£alse_nearest from the TISEAN package. However, because 
our data are both short and noisy, we have to skip the second 
criterion, which omits points being farther than the standard de¬ 
viation of the data divided by the threshold. Otherwise we would 
have very small number of points to consider and the results are 
meaningless. Therefore we are forced to include also points, 
which are not very close neighbours (even when they are the 
nearest ones), because due to the noise and the limited amount 
of data, there is not enough really close points. We also have to 
be careful about the size of the threshold Rtoi. which cannot be 
very high, otherwise almost no nearest neighbour distance can 
grow by this factor due to a finite size of the attractor. This can 
slightly affect the results. 

We choose the threshold Rtoi = 10. The corresponding 
false nearest neighbour ratio for the set of IGR’s observations is 
given in Fig. |A.2| We set the embedding dimension m = 10, in 
which for almost all observations the false nearest neighbour ra¬ 
tio drops to zero. However, even m - 6 could be a good choice, 
because for this value there is only 1% or less of false neigh¬ 
bours. 


These two methods are used for guessing some appropriate 
values of the two parameters, but the results of the recurrence 
analysis do not strongly depend on it. This is supported by|Thiel| 


et al. (|2004)l, who claimed that some of the dynamical invariants 


estimated by recurrence analysis including second order Renyi 
entropy K 2 do not depend on those parameters. 
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Fig. A.3. RP for several observations of IGR J17091-3624 from Table 1. row: left - SIMS/HIMS state 02-00 (left upper comer in blue) and 
p state 06-02 (right lower comer in red). Long diagonal lines appear in the red RR Right - The zoom of preceding plot. The diagonal lines are 
affected by the noise in such a way, that the apparently long lines are composed of shorter square-shaped pieces.Next rows: In bottom red halves 
of each plot RP is plotted for observations 04-01 (2. row left, IVS state), 04-02 (2. row right, p-state), 05-02 (3. row left, SIMS state) and 07-01 
(3. row right, p-state). In green colour in the left upper half of each plot one of the corresponding surrogates is shown. 
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Having the values for embedding delay and dimension de¬ 
termined we can now proceed with the recurrence analysis. The 
basic object of the analysis is the recurrence matrix defined as 
follows: 


R,;;(e) = ©(e- II X,- - X,. ||), i, j ^ I, (A.2) 

where x, = x(f,) are (N) points of the reconstructed phase tra¬ 
jectory and 0 is the Heaviside step function. The matrix thus 
contains only I’s and O’s and can be easily visualised by plot¬ 
ting a dot at the coordinates i, j whenever R,j(e) = 1, which 
is called the recurrence plot (RP). From the definition ( |A.2| l it is 
obvious, that the RP is symmetrical with respect to the main di¬ 
agonal, which is formed by dots (for i - j) and which is omitted 
in further computations. 

The comparison of the RPs for two observations, 02-00 
(SIMS/HIMS - upper left corner on the left panel in first row) 
and 06-02 (p-class - lower right corner of the same panel) is 
given in Fig. |A.3| Because the plot is symmetrical with respect 
to the main diagonal, we show only one half of the RP for each 
observation in the upper/bottom trilateral half of the plot. Differ¬ 
ent geometrical structures contained in the RP could be seen in 
these two examples. The structures can be quantified based on 
the statistical properties of the recurrence matrix. 

The most important entity for the recurrence analysis is the 
existence of long diagonal lines, which indicate periodic (and 
regular) behaviour. Such structures are much more prominent 
in the RPs of p observations. The quantification of the visual 
information is contained in the histogram of diagonal lines of a 
certain prescribed length I, 


N Z -1 

P{e,l) = (A.3) 

ij=i 


k=0 


Based on the histogram more other quantifiers can be defined 
(Marwan et al. 2007 1 . Here we will focus mainly on Lmax = 
max({f,j" j), which is the length of the longest diagonal line (ex¬ 
cept of the main diagonal), and on the estimate of the second- 
order Renyi’s entropy K 2 . This entropy, also called correlation 
entropy or correlation exponent, collision entropy or just Renyi’s 
entropy ( Grassberger||1983 i, is defined for the phase trajectory 
x(f) as 


K2 - - lim lim lim — In 

A/—>0 e —>0 A/—>00 l^t ^^ 


P\,. 




(A.4) 


where the phase space is coarse-grained into m-dimensional 
hyper-cubes of size e. Here, describes the probability that 

the first point of the trajectory x(f = At) is inside the box i\, the 
second point x(f = 2Af) is in the box i 2 and so forth to the /-th 
point being in the box //. It can be shown ( |Marwan et al.|2007 1 , 
that the square of this probability is connected with the existence 
of a diagonal line in the RP (the diagonal line means, that the 
trajectory goes twice through the same sequence of phase-space 
boxes). 

More precisely, K 2 is related with the cumulative histogram 
of diagonal lines pde, 1), describing the probability of finding a 
line of minimal length I in the RP, by the relation 

Pc{e,r) ~ (A.5) 


hence we can estimate the value of K 2 as the slope of the loga¬ 
rithm of the cumulative histogram versus / for constant e (see an 
example in Fig. A.4 1 . This estimation holds for the limit I —> 00 , 
particularly for I > At. The important property of K 2 entropy 



Fig. A.5. Renyi’s entropy K 2 versus recurrence rate for the obser¬ 
vation 06-02 for three different embedding dimensions marked in the 
figure. 



Fig. A.6. Dependence of the longest diagonal line given in points 
on the recurrence threshold e for the observation 05-02 (thick red line). 
The ensemble of 100 surrogates are shown in grey and the ensemble of 
shuffled surrogates are shown in cyan. 


is that it is the lower estimate of the sum of positive Lyapunov 
exponents of the system, hence it is positive for chaotic systems 
( |Marwan et al.|2007l l. 

Considering formula ( |A.5| l for two different values of thresh¬ 
old 62 -^ 1 + Ae, we can also estimate the value of the correlation 
dimension D 2 by 


D 2 (e) = In 


l pde\,l) 

KPcieiJ) 



(A.6) 


However, for obtaining reasonable results for D 2 much more 
data points are needed than for estimating K 2 ( |Marwan et al.| 
2007| l. Because we are limited by the length of individual ob¬ 
servations, which are not longer than several thousand seconds 
yielding twice data points for our chosen binning time df = 0.5 s, 
we were unable to determine the value of D 2 by means of the 
recurrence analysis. 

Besides the short duration of the observation, another diffi¬ 
culty is the presence of noise, which affects the structures con¬ 
tained in RP. This is documented by the right panel in first 
row of Fig. |A.3| which shows the zoom of the left panel for 
t e (50,200) s. Here the apparently long lines in fact consist 
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Fig. B.l. The same as in Fig. |A.6| for the observation 05-01. The 
condition I > At on lines’ length used for A' 2 -estimation is indicated by 
black horizontal line. 


of shorter and broader pieces separated by gaps. The gaps ap¬ 
pear, when the noise takes the “should-be recurrence point” out 
from the e-neighbourhood, and on the other hand, the edges on 
the lines appears, when noise brings the “should-not-be recur¬ 
rence point” inside. Hence, there is more shorter diagonal lines 
and less longer ones, thus the slope of the histogram is steeper. 

The most usual way to deal with the noise in data is to take 
the radius of the neighbourhood big enough, so that it covers 
the noisy perturbations of the position, which according to [Thiel 
|et al.j ( |2002| l; [Marwan et al.j ( |2007| l is e = 5cr, where cr is the 
standard deviation of the noise. However, a part of the noise 
is caused by the measurement process and another part stems 
from the stochastic processes behind the emission of the radia¬ 
tion. Moreover, the amount of noise and its standard deviation 
is unknown, but it is quite large comparing to the standard de¬ 
viation of the whole data set (see Fig. [^l. Hence, such a choice 
of threshold would not be appropriate, because it would be huge 
and almost every point would become the recurrence point. 

We consider such thresholds, which yields the recurrence 
rate (RR - ratio of the recurrence points to all points of the recur¬ 
rence matrix) within 1% to 25% and we have to keep in mind, 
that due to the noise, the slope of the cumulative histogram and 
K 2 is overestimated. An example of the decreasing dependence 
of the /r 2 -estimate on the recurrence rate with At -7 s, and three 
different embedding dimensions mi = 10, m 2 = 20 and m 3 = 30, 
is given in Fig. |A.5| 


Appendix B: Comparison between the surrogate 
data and observations of iGR J17091 

We first study the appearance of long diagonal lines in RP. For 
the real data set and its surrogates we compare the dependence 
of Tmax on chosen e. We obtain three different types of output. 

For observations from states other than p and for the used 
range of thresholds we do not see any lines longer than At, hence 
it is impossible to estimate K 2 . Moreover, in Fig. A .6 we can see 


that the length of the longest line for the data (thick red line) and 
its surrogates (grey lines) do not differ. We conclude that the 
origin of these observations is sufficiently well described by the 
assumption of linearly autocorrelated Gaussian noise. We can 
also try to soften our null hypothesis and use the shuffled surro¬ 
gates. It turns out, that all the observationswith short lines only 
are also perfectly consistent with this null hypothesis, thus here 
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Fig. B.3. Comparison of recurrence rate for observation 06-02 (value 
is indicated as the horizontal thick red line) and its hundred surrogates 
(number of surrogate is on the x-axis) computed for the same recurrence 
parameters (e = 2.8, ra = 10, At = 14dr). 


is no evidence of any dynamics at all (in Fig. |A. 6 | the shuffled 
surrogates are plotted by cyan color). 

Next group of observations are those, which do contain 
longer diagonal lines (usually for higher e), but the values of 
Lmax for data and lAAFT surrogates do not differ notably. Within 
our sample of observations of IGR J17091-3624 only the obser¬ 
vations 05-01 and 06-03 can be assigned to this subgroup, the 
example for 05-01 can be seen in Fig. EH However, compar¬ 
ing the data with the shuffled surrogates (cyan lines in Fig. ED’ 
we find a striking difference, meaning that those observations 
cannot be described as temporally independent identically dis¬ 
tributed random data. Hence in this case our analysis suggests 
that the data comes from a linear dynamical process (linearly 
autocorrelated Gaussian noise). 

The rest of p observations show stronger difference between 
the data and its surrogates, is higher than for all sur¬ 
rogates for wide range of e. This behaviour indicates non-linear 
features of the background dynamics. 

In order to quantify these results, we compute the estimate of 
K 2 for observations with enough long lines and we will use this 
quantity as our discriminating quantity. The computation goes 
as follows. 

For a chosen value of e and the chosen embedding dimen¬ 
sion and delay {m — 10 and At -7 s) we compute the recurrence 
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Fig. B.4. Comparison of estimates of K2 for observation 06-02 (indi¬ 
cated by horizontal thick red line) and its hundred surrogates (number of 
surrogate is on the z-axis) computed for the same recurrence parameters 
(e = 3.25,m= 10,Af= 14d?). 


matrix and the cumulative histogram for the data and for its sur¬ 
rogates. Usually the recurrence rate for the surrogates is lower, 
which means that the shuffled order of points yields more dis¬ 
tant points after the reconstniction of the phase space trajectory 
(see Fig. |B.3[ ), however for higher recurrence thresholds this dif¬ 
ference is smeared. Afterwards, we estimate K 2 from the slope 
of the cumulative histogram for the central part of the histogram 
which satisfies / > Af, Pc(e,l) > Nmin and I < O.lNdt. We also 
require that there is at least 5 points of the cumulative histogram, 
which satisfy our criteria and could be used for the linear regres¬ 
sion. The value of Amin has to be chosen with respect to the 
length of the time series and the total number of lines. We usu¬ 
ally choose Amin = 100. The example of cumulative histogram 
computed for the observation 06-02 with e - 2.8 is given in 
Fig. |A.4[ where the used region of the histogram is marked by 
two vertical dashed lines. 

We execute the same analysis with the same parameters on 
the set of surrogates and we obtain set of values (see 
Fig. |B.4[ ). According to [Theiler et al.| ( [1992) 1 we define a sig- 
nihcance of the obtained result as a weighted difference between 
the estimate of Renyi’s entropy for the observed data and 
the mean of estimates of Renyi’s entropy for ensemble s of the 

the 


A.5 


surrogates However, as can be seen from relation 
quantity is non-negative and behaves in exponential manner 
(for 10-times longer diagonal lines, K 2 drops down by one or¬ 
der). Therefore we first introduce the quantity Q - ln(A2), which 
scales linearly, and then dehne the significance stemming from 
the estimate of K 2 as 


SKiie) - 


cr 


(B.l) 


where (jQiar is the standard deviation of the set (see 

Fig.|B;g. 

When the observed data yields but the RP of its surro¬ 
gates does not contain enough long lines, it is impossible to com¬ 
pute . Obviously this should contribute to the significance 
in the case, that is lower than A™ . However, the values of 
gsuiT jjQj available. We set the significance of the case, when 
all A™‘^ surrogates have only short lines, to be .Ssi = 3. 

On the other hand, if is higher than A™" , then some of 
the surrogates have lines longer than the observed data. With that 
thecoincident existence of surrogates, which do not have long 


Fig. B.5. Comparison of Q for observation 06-02 (indicated by 
horizontal thick red line) and its hundred surrogates (number of sur¬ 
rogate is on the jc-axis) computed for the same recurrence parameters 
(e = 3.25, m = 10, A? = 14df). The mean value Q and the standard 
deviation of the ensemble of surrogates is indicated by horizontal lines. 


lines at all, lessens the signihcance. This fact is incorporated in 
our procedure in the way that we add or subtract the two sig¬ 
nificances <Ssi and Sk 2 according to the sign of the difference 
e°bs(e) - 2--(e). 

Finally, we obtain the total significance as the weighted dif¬ 
ference of signihcances of these two cases, namely 

■5(') = - sign(a*(e) - (B.2) 

where Aji is the number of surrogates, which have only short 
lines, and A^“" = Nsg. +Asi is the total number of surrogates. Our 
definition of signihcance yields negative values for some cases, 
e.g. when and Asi = 0. Since we expect, that the 

real non-linear dynamical system should contain more long lines 
than the artihcial surrogates, only positive values of signihcance 
serve as the indication of non-linear dynamics. In fact, we use 
the value S - 1.5 as the limit, above which we consider the 
observation to be produced by non-linear dynamical system. 

The last ambiguous step in the calculation of the signihcance 
is the choice of the recurrence threshold e. There are two pos¬ 
sible ways how to treat this issue. One way is to choose some 
recurrence rate for the observation (e.g. RR ~ 10%), choose the 
threshold appropriately and compute the signihcance. This ap¬ 
proach is suitable for quick analysis, however it is not very accu¬ 
rate. In this paper we adopted the other way, which lies in hnding 
the range of e, for which RR 6 (1%, 25%), computing the signif¬ 
icance for a linearly spaced set of e in this interval (typically, we 
have used ~ 10-40 values of e in dependence on the number 
of points of the data set) and establishing the hnal decisive value 
as the average of the obtained signihcances S - 
If for some e the observed data set does not have enough long 
lines for the K 2 estimation, this value is not taken into account in 
the evaluation. 

In Table [^ we summarized the length of the longest diagonal 
line (expressed as the number of points) for three different recur¬ 
rence rates_ RR~5%, 10% and 20%, the resulting averaged sig¬ 
nihcance S and the number of values of e used for the averaging. 
All p observations except of 05-01 and 06-03 (recognized earlier 
as special group with respect to Lmax behaviour) have S > 1.5. 

In the same way we can also compute the signihcance .Sshf 
for the other null hypothesis using the set of shuffled surro¬ 
gates. This quantity serves as the indicator of departure from 
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stochastic behaviour. Typically, we compute this quantity only 
for those observations, which have enough long lines but show 
non-significant result with the lAAFT surrogates. This is be¬ 
cause the data sets which have significant result for the non¬ 
linear dynamics should automatically have significant result for 
non-stochastic nature. However, for the confirmation of this as¬ 
sumption in Section]^ we also provide several examples of this 
quantity for observations displaying the non-linear features. 


Appendix C: Testing the method with simuiated 
time series 

In general, our method can be applied to different kinds of time 
series, which are the product of any dynamical system, not only 
to X-ray lightcurves. For testing it is desirable to apply the 
method on time series, whose nature is known and which has 
been studied by other means. Because the dynamical system be¬ 
hind the observed lightcurves is a priori not known, we rather 
use numerical data. We chose the numerical time series, which 
describe the motion of geodesic test particle in the field of static 
black hole surrounded by a massive thin disc. 

The details about the testing procedure will be given in sep¬ 
arate paper ( Sukova & Janiuk]|2015, in preparation| l. Here we 
only summarise the main results. 

Our analysis performed on the numerical time series revealed 
what the recurrence analysis provides for regular and chaotic 
trajectories of the complicated non-linear system, whose phase 
space consists of regular KAM tori nested in the chaotic layers. 
We have seen that even though governed by the same global set 
of equations, the evolution of the orbits in the regular and chaotic 
parts of the phase space can differ significantly, which is also re¬ 
flected by the patterns in the recurrence plots. In agreement with 
the fact, that regular orbits in conservative systems have zero 
Lyapunov exponents, the estimate of K 2 is much lower for the 
regular trajectory than for its chaotic counterpart. 

Despite the expectation, that regular motion should not ex¬ 
hibit non-linear behaviour in our analysis, the significance for 
the regular orbit is quite high, however still almost twice lower 
than for the chaotic one. We argue that the reason for this is the 
way how the surrogates data are constructed, which means that 
they have exactly the same value distribution but they reproduce 
the spectrum only approximately, depending also on the avail¬ 
able length of the data set. In case of the regular motion, very 
narrow peaks are in the spectrum and the error in reproducing 
such spectrum causes the very long diagonal lines to be broken. 

However for increasing strength of the added noise, the sig¬ 
nificance of the regular orbit drops down very quickly, while the 
chaotic orbit shows solid signs of non-linear behaviour even in 
the case, when the added noise has the same variance as the orig¬ 
inal data. Yet, even for high noise levels, L^ax for the regular 
orbit is higher or similar as for the chaotic orbit. Hence, obser¬ 
vations which show non-significant result with the lAAFT sur¬ 
rogates, but very significant result for shuffled surrogates, and 
contain long diagonal lines compared to the time delay Af, can 
possibly be explained as regular motion. Another possibility is, 
that such observations cover the part of the chaotic trajectory, 
which is called “sticky orbit”. Because the observations cover 
only limited time interval, it is not possible to distinguish be¬ 
tween these two cases. 
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